Testing Proposals for the Yang-Mills Vacuum Wavefunctional 
by Measurement of the Vacuum 
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We review a method, suggested many years ago, to numerically measure the relative amplitudes of the true 
Yang-Mills vacuum wavefunctional in a finite set of lattice-regulated field configurations. The technique is 
applied in 2+1 dimensions to sets of abelian plane wave configurations of varying amplitude and wavelength, 
and sets of non-abelian constant configurations. The results are compared to the predictions of several proposed 
versions of the Yang-Mills vacuum wavefunctional that have appeared in the literature. These include (i) a 
suggestion in temporal gauge due to Greensite and Olejnflc; (ii) the "new variables" wavefunction put forward 
by Karabali, Kim, and Nair; (iii) a hybrid proposal combining features of the temporal gauge and new variables 
wavefunctionals; and (iv) Coulomb gauge wavefunctionals developed by Reinhardt and co-workers, and by 
Szczepaniak and co-workers. We find that wavefunctionals which simplify to a "dimensional reduction" form 
at large scales, i.e. which have the form of a probability distribution for two-dimensional lattice gauge theory, 
when evaluated on long-wavelength configurations, have the optimal agreement with the data. 
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I. INTRODUCTION 

Most of the key non-perturbative properties of non-abelian 
gauge theories, such as the static quark potential, the chi- 
ral condensate, and the topological charge density, are actu- 
ally properties of the vacuum of the quantized theory. In the 
Hamiltonian formulation, the vacuum state is the ground state 
wavefunctional of the Hamiltonian operator, and all of the ex- 
cited states of the theory, i.e. the mesons, baryons, and, in 
a pure gauge theory, the glueballs, are simply small excita- 
tions on top of that underlying ground state. For this reason, 
knowledge of the Hamiltonian ground state wavefunctional 
could be essential in understanding the infrared properties of 
a non-abelian gauge theory. 

Proposals for the ground state of pure Yang-Mills theory go 
back over thirty years [1,2]. However, with only a few excep- 
tions [3-7], very little work was done in this area after those 
initial efforts. In recent years, however, there has been a mod- 
est revival of interest in this area, and a number of plausible 
suggestions for the vacuum state have been advanced. These 
proposals will be described, along with their motivations, in 
the next section. Briefly, there are suggestions which have 
been put forward in temporal gauge [8], in Coulomb gauge 
[9-13] and, in 2+1 dimensions, in terms of gauge-invariant 
"new variables" [14]. Since these suggestions differ in vari- 
ous ways, it would be interesting to know which (if any) is the 
true vacuum state, or at least a reasonable approximation to 
the true vacuum state. 

In this article we will apply an old method [15-17] for mea- 
suring, via lattice Monte Carlo simulations, the relative mag- 
nitudes of the true Yang-Mills wavefunctional in any given set 
of lattice gauge field configurations. The evaluations will be 
carried out for two types of lattice configurations: non-abelian 
constant gauge fields of varying amplitudes, which are con- 
stant in space but noncommutative [Ui,Uj\ 7^ 0, and abelian 
plane waves of various amplitudes and wavelengths, which 



are abelian in the sense that [Ui,Uj\ = 0. The results are com- 
pared to the corresponding values obtained in each of the pro- 
posed vacuum wavefunctionals. The method can be applied in 
any number of space-time dimensions, but here we will work 
exclusively in 2+1 dimensions, since the new variables pro- 
posal [14] is formulated only in that case. 

In section II below we will introduce and motivate each 
of the wavefunctionals to be tested. Section III reviews the 
method for measuring the true vacuum wavefunctional, and 
section IV compares the results obtained by this method with 
the predictions of each of the proposed ground states. Our 
conclusions are in section V, and some numerical details are 
found in the appendix. 



H. VACUUM STATE PROPOSALS 

The Yang-Mills Hamiltonian operator takes on its simplest 
form in temporal gauge, namely 
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in the continuum theory in D + 1 dimensions, and 
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on the lattice, where the sums are over links I and spatial pla- 
quettes p, respectively. Physical states in temporal gauge must 
obey the Gauss law constraint Df'EpV = 0, or more explicitly 



(3) 



which implies that physical states must be invariant under in- 
finitesimal gauge transformations. The Gauss law constraint 
in temporal gauge is a mixed blessing in the search for an 
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approximate ground state. On the one hand, gauge invariance 
can be seen as an aid in selecting a good ansatz for the vacuum 
state. On the other hand, by severely limiting the choice, cer- 
tain states which are perfectly acceptable in Coulomb gauge, 
and which may be much more amenable to an analytical treat- 
ment, must be discarded in temporal gauge. A very impor- 
tant relation, for our purposes, is the equality of the vacuum 
wavefunctionals in temporal and Coulomb gauge (see, e.g., 
ref. [18]), 



Coul [A}=^ mp [A] 
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when evaluated on gauge fields satisfying the Coulomb gauge 
condition V • A = 0, and which also lie in the first Gribov re- 
gion. Since our numerical method, to be described in the next 
section, will generate the relative amplitudes of vacuum wave- 
functionals in temporal gauge, in any finite set of gauge field 
configurations, we will be able to check proposals in Coulomb 
gauge by ensuring that the given set satisfies the Coulomb 
gauge condition, and lies within the first Gribov horizon. 

The ground state wavefunctional is known in two limits: 
the free-field g 2 = limit, and also at strong lattice couplings 
g 2 3> 1 ■ In the free-field limit, in either Coulomb or temporal 
gauge, 



«P [A]=exp 
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while in the strong-coupling limit, in SU (N) gauge theory, it 
has been shown that [19] 
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to leading order in l/g 2 . It was suggested long ago in ref. 
[1], by one of the present authors, that the Yang-Mills vacuum 
wavefunctional in 3+1 dimensions might have the form 
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when evaluated on sufficiently long-wavelength, slowly vary- 
ing field configurations. This wavefunctional has the property 
of dimensional reduction: If we write 
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then R [A] has the form of the Euclidean Yang-Mills action in 
one lower dimension (three dimensions, in this case). It is 
clear that the strong-coupling vacuum state (6) does, in fact, 
have this property. 

The dimensional reduction vacuum (7) in 3+1 dimensions 
is confining, i.e. 



W(C) = 0P o |Tr[E/(C)]|%) 
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if and only if Yang-Mills theory in three Euclidean dimen- 
sions has that property, where U(C) is a Wilson loop holon- 
omy around the planar, spacelike loop C. Of course we have 



good reasons to believe that Yang-Mills theory is confining 
in three Euclidean dimensions. It was noted by Halpern [2] 
that a dimensional -reduction vacuum state in 2+1 dimensions 
must be confining, since Yang-Mills theory in two Euclidean 
dimensions is known to confine. Dimensional reduction was 
also suggested somewhat later, on rather different grounds, 
by Ambjorn, Olesen, and Peterson [20, 21]. These authors 
were the first to make the connection between dimensional re- 
duction and the property that has come to be known [22] as 
Casimir scaling. Strong evidence for Casimir scaling at inter- 
mediate distance scales was found in [23]. 

On the other hand, the dimensional reduction wavefunc- 
tional cannot be correct as it stands, because the short-distance 
structure is completely wrong. For example, equal-time two- 
point correlators in D + 1 dimensions, at short distances, 
cannot be identical to short-distance two-point correlators 
in D Euclidean dimensions; the singularity structure in the 
approach to zero separation would be wrong. In general 
one would expect that the vacuum state evaluated on short 
wavelength configurations would agree with the perturbative 
ground state, whose zeroth order approximation is given by 
(5). 

There are other reasons, apart from short-distance singular- 
ity structure, that dimensional reduction cannot be exact even 
for infrared physics. Dimensional reduction from 2+1 to two 
Euclidean dimensions would imply a non-vanishing string 
tension, and perfect Casimir scaling, for any color group rep- 
resentation. This cannot be right in 2+1 dimensions, because 
of color screening. 1 As argued in ref. [8], it is quite plausi- 
ble that color screening is achieved by small corrections to the 
dimensional reduction form. 

Another argument against exact dimensional reduction 
from 3+1 to three Euclidean dimensions was raised in refs. 
[26, 27], which pointed out that this reduction would imply a 
match between the equal-time Coulomb gauge gluon propa- 
gator in 3 + 1 dimensions, and the Landau gauge propagator 
in D = 3 Euclidean dimensions. It was shown in the same 
references that these propagators actually do agree quite well 
in a range of low and intermediate momenta around 1 GeV (a 
range which is relevant for phenomenology), but the equiva- 
lence cannot hold in the far infrared. 

For all of these reasons, a purely dimensional reduction vac- 
uum wavefunctional is clearly inadequate. Corrections are es- 
sential, and what is really required is an approximation to the 
vacuum state which holds at all distance scales. There are now 
a number of proposals, which may or may not obtain the di- 
mensional reduction form in some limit, but which do claim 
to approximate the ground state at all length scales. These we 
will briefly review. 



1 For this reason it is useful to consider ^-string tensions, associated with 
quarks in completely antisymmetric representations, whose color charge 
cannot be screened to a lower dimensional representation by gluons. The 
current evidence [24] in 2+1 dimensions is that the leading corrections to 
the JV = °° result are of order 1 /N , as in Casimir scaling, rather than 1 /N 2 , 
as in the competing Sine Law proposal. For a recent discussion of ir-string 
tensions in the context of the large-W expansion, cf. [25]. 
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A. Temporal gauge 

It was suggested in ref. [8] that the Yang-Mills ground state 
wavefunctional, in D = 2 + 1 dimensions and in temporal 
gauge, is approximated by 2 
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where B" = F" 2 , D 2 is the covariant Laplacian, Ao is the low- 
est eigenvalue of — D 2 , and m 2 is a parameter which vanishes 
as g — > 0. The motivation was to find the simplest possible 
gauge-invariant expression which would agree with the free- 
field (5) and dimensional reduction (7) wavefunctionals in the 
appropriate limits. In support of this conjecture, it was found 
that ^ao 

1. solves the Yang-Mills Schrodinger equation in the 
strong-field, zero-mode limit; 

2. confines if the mass parameter m > 0, and that m > 
seems to be energetically preferred; 

3. produces results for the mass gap, the Coulomb gauge 
ghost propagator, and the color Coulomb potential, 
which are in rather good agreement with results derived 
from standard lattice Monte Carlo simulations. 

The subtraction of Ao is essential, and was introduced be- 
cause —D 2 has a positive semi-definite spectrum, and in gen- 
eral the lowest eigenvalue tends to infinity for typical vacuum 
configurations in the continuum limit. This fact is obvious 
perturbatively, and is confirmed numerically. Without the sub- 
traction (and this was the form originally suggested by Samuel 
[6]), the kernel joining B a (x) and B b (y) in (10) effectively van- 
ishes in the continuum limit, and the corresponding string ten- 
sion would be infinite. In contrast, the spectrum of — D 2 — Ao 
is well-behaved, and not far from that of the free-field Lapla- 
cian operator — V 2 [8]. 

If one drops all components of the vector potential apart 
from the zero mode (analogous to the "minisuperspace" ap- 
proximation in quantum cosmology), then the Lagrangian and 
the Hamiltonian operators are simply 



d,A k ■ d t A k - (Aj x A 2 ) ■ (Ai x A 2 ) 
d,A k ■ d,A k - (Ai x A 2 ) • (Ai x A 2 ) 



L =^ 2 J d2x 
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H=-4=^— + ^- 1 (A l xA 2 )-(A l xA 2 ) 7 (11) 
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2 A factor of g has been absorbed into the definition of the gauge field, so 
that Ak has units of inverse length. This accounts for the overall factor of 
1/g 2 in the exponent of the wavefunction. 



where V is the volume of 2-space, and the cross-product and 
dot-product are defined with respect to SU (2) color indices. 
Solving for the ground state is a problem in quantum mechan- 
ics, rather than quantum field theory, and to leading order in 
1 /V the solution is 



= exp 
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Now in the region of parameter space where the zero mode is 
much larger than all other modes, the covariant Laplacian is 
approximated by 
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and m 2 is negligible. It is then found, after some algebra, that 
the proposed wavefunctional (10) reduces to the zero-mode 
solution (12). 

Dimensional reduction follows by expanding the B-field in 
eigenmodes <j>" of — D 2 . Then the part of the wavefunctional 
that depends only on the low-lying modes, with eigenvalues 
A„ — Ao <C m 2 has the form of the dimensional reduction wave- 
functional (7), with 11 = 1/ m. If we assume that the asymp- 
totic string tension is due to the low-lying modes, then cal- 
culation of the string tension is simply an exercise in two- 
dimensional Yang-Mills theory, and the result is 
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If we turn this around, and write m = 16<7/(3g 2 ), then we have 
a complete proposal for the vacuum wavefunctional, although 
the string tension must be supplied as an input. 

A method for obtaining equal time expectation values 
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by numerical simulation, with a suitable lattice regulariza- 
tion, was also introduced in [8], and applied to calculate the 
mass gap. The Coulomb gauge ghost propagator and color 
Coulomb potential were derived via numerical simulation of 
in [28], by the method of generating thermalized lattice 
configurations from the ^Vqq distribution, and then transform- 
ing these configurations to Coulomb gauge. The results, as 
already mentioned, were in very good agreement with those 
obtained from standard lattice Monte Carlo simulations. For 
details, we refer the reader to the cited references. 



B. New variables 

While the temporal gauge ground state can be credited with 
some numerical success, it remains an educated guess, and 
requires the string tension as an input. A more ambitious pro- 
gram in 2+1 dimensions, which aims to calculate both the 
Yang-Mills vacuum state and the string tension analytically, 
was initiated by Karabali, Kim, and Nair [14], and has been 
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further developed by Karabali and Nair in a series of papers, 
cf. [29] and references therein. 

The starting point in the Karabali, Kim, Nair (KKN) ap- 
proach is temporal Aq = gauge, and the remaining two 
components of the A-field are combined into a complex field 
A — (A\+ iA^i)j2, related to a matrix-valued field M via 



A = -(<?-M)Ar 1 , A = M t ~ 1 o*jM t , 



(16) 



where z = x\ — ixi, and z = x\ + ixz are the usual holomorphic 
variables in the complex plane. The matrix-valued field M 
takes values in the group SL(2,C), and transforms covariantly, 
M — > GM, under a gauge transformation G. This field can be 
used to define gauge-invariant field variables 



ye = M^M 

a _c A d.ye 



ye- 



(17) 



where Ca is the quadratic Casimir in the adjoint representa- 
tion. In terms of these gauge invariant variables, the Hamilto- 
nian becomes 

H KKN = T + V, (18) 

where T is derived from the E 2 term in the standard Hamilto- 
nian 

8 



T = m 
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and also 
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Inner products are evaluated with respect to the integration 
measure 

(«Pi|«p 2 ) = [ dn(ye) e 2CASwzw ^ )x ¥* 1 (ye)^ 2 (ye) ,(23) 



where djj. (ye) is the Haar measure, and Swzw is the Wess- 
Zumino-Witten action. 

Although the new field variable J? is gauge invariant, the 
Hamiltonian Hkkn is invariant under local holomorphic trans- 
formations h(z), under which J! transforms like a connection 



-*h fhT l + —dhhr l 



(24) 



and all physical states in the new variables approach, 

must be invariant under this local transformation. In this 



sense, the new variables approach trades the local gauge in- 
variance constraint (the Gauss law) in temporal gauge for in- 
variance under local holomorphic transformations. 

Expressing the ground state as ^oij?] = ex Pi~R[y?))> 
KKN find an expression for R\ c f\ which is bilinear in ^ , 
namely 



Vkkn = ^exp 
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where the second line is the new variables state converted back 
to usual variables. KKN assume that the dimensional reduc- 
tion form is obtained for long-wavelength configurations by 
simply dropping —V 2 in the kernel, i.e. 



KKN 



— > ^Kexp 
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and then the string tension for a spacelike Wilson loop is ob- 
tained from solving Yang-Mills theory in two Euclidean di- 
mensions, with the result 



(27) 



Very remarkably, this value is within a few percent of the value 
found by Bringoltz and Teper [30] in lattice Monte Carlo sim- 
ulations of the 2+1 dimensional theory, after careful extrapo- 
lation to the continuum limit. 3 



C. A hybrid wavefunctional 

The problem with ^kkn is that, in terms of new variables, 
it is not holomorphic invariant, and in terms of the usual vari- 
ables (second line of (25)) it is not gauge invariant. Therefore 
^kkn, as it stands, is not a physical state. Of course, KKN 
do not claim that x ¥kkn\J ! \ in eq. (25) is exact, and pre- 
sumably gauge and holomorphic invariance requires consid- 
eration of contributions to R\J?\ involving higher powers of 
J" . However, ignorance of the gauge/holomorphic-invariant 
wavefunctional calls into question the assumed dimensional 
reduction form (26), which was required for the successful 
prediction of the string tension. For example, suppose we 



3 Recently some corrections to a have been calculated [29], and they are 
quite small. At present it is not entirely clear why the correction is so 
small, since there is no obvious small expansion parameter in this approach, 
and the corrections involve a sum of rather large (positive and negative) 
contributing terms, which for some reason nearly cancel. 
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assume that higher powers of @ in the expansion of R\J?\ 
would have, as its main effect, the conversion of the ordinary 
Laplacian into a covariant Laplacian; i.e. in the usual variables 



= c/^exp 
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(28) 



In that case, for configurations which are non-abelian 
([A;t,A v ] 7^ 0) in general, dropping — D 2 is invalid even for 
configurations which vary very slowly compared to the length 
scale 1 /g 2 , and indeed is invalid even for configurations which 
have no spatial variation whatever. As we have remarked 
above, in connection with ^GO; the covariant operator — D 2 
has a positive semi-definite spectrum, and for typical lattice 
configurations the lowest eigenvalue diverges in the contin- 
uum limit. In that case, rather than replacing — D 2 by zero to 
obtain the dimensional reduction result, one should replace it 
by infinity! This is obviously nonsense. 

Assuming that the KKN wavefunctional applies to abelian 
configurations ([A r ,A v ] = 0), the corresponding vacuum state 
for more general configurations is still a mystery; one can 
only guess what the gauge and holomorphic invariant com- 
pletion of ^kkn might be. But the gauge-invariant comple- 
tion is essential, if one is going to invoke dimensional reduc- 
tion to compute the string tension. At this stage there are an 
infinite number of possibilities, and the validity of the KKN 
prediction for the string tension depends on which of these 
possibilities is the correct one. One possible approach is to 
retain ^kkn for abelian configurations, and ask for the sim- 
plest gauge-invariant generalization which would lead to the 
dimensional reduction form (26). Then it is natural to merge 
features of ^go an d ^kkn into a conjectured "hybrid" form 
for the ground-state wavefunctional 



hybrid 
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which we will include in our numerical tests below. 

An alternative approach has been followed by Leigh, Minic, 
and Yelnikov (LMY) [31], who begin with the ansatz 



Wlmy = exp 



d 2 xd 2 ydf a (x)K xy (L)djf a (y) 



(30) 

whereL= —A/m 2 , andAis the holomorphic -co variant Lapla- 
cian. They then derive and solve a differential equation for 
K{L), where L is treated as a number, rather than an operator, 
and by solving this equation they arrive at 



K(L) 



1 J 2 (4-Vl) 



(31) 



where J\ 2 are Bessel functions. By construction, the LMY 
proposal is a physical state. If the infrared limit means L — >• 0, 



then K—¥\, and has the dimensional reduction form (26), 
leading to the same prediction for the string tension. Leigh 
et al. also obtain predictions for the glueball mass spectrum 
in 2+1 dimensions, which appear to be in good agreement 
with standard lattice Monte Carlo results. The reservation in 
this case is that the LMY approach assumes a certain operator 
identity (eq. (56) of ref. [31]) whose validity, in our opinion, 
is questionable. It would nevertheless be interesting to test 
^lmy numerically, but unfortunately it is not clear to us that 
the method we will use in this article could be easily applied 
to the LMY proposal. 

D. Coulomb gauge 

In Coulomb gauge, after resolving Gauss' law, eq. (3), one 
obtains the Yang-Mills Hamiltonian [32] in terms of the trans- 
verse components of the gluon field, V • A = 0, 

- 1 [A] U1J [A] + B?B?) +H C (32) 



H 



d D x{ 



H C = S -J d D xd D yf-'[A]p a (x)f[A]F« b (x,y,[A])p b {y), 

where 11" (x) = 8/i8A"(x) is the canonical momentum (elec- 
tric field) operator and 

/A]=Det(-D-V) (33) 

is the Faddeev-Popov (FP) determinant (this should not be 
confused with the variable J" (x) in the KKN approach). Fur- 
thermore 

p a (x) 

is the color charge of the gluons and 



e abc A b n c 



(34) 



F ab {x,y,[A\) 



-D-V)- 1 (-V 2 )(-D-V)- 



x.a\y,b 



(35) 



is the so-called Coulomb kernel. The gauge fixed Hamilto- 
nian eq. (32) is highly non-local due to the Coulomb kernel, 
eq. (35), and due to the FP determinant, eq. (33). In addition, 
the latter occurs also in the functional integration measure of 
the scalar product of Coulomb gauge wavefunctionals 



(36) 



{Wx\OWi)= DAf[A] W l[A]Oy 2 [A\. 



Any normalizable state, expressed as a functional of the trans- 
verse gauge field, is a physical state in Coulomb gauge. This 
means in particular that a wavefunctional which is Gaussian in 
the gauge field may be a viable proposal for the ground state. 
Unlike the GO and KKN/hybrid proposals, such a state cannot 
have the dimensional reduction property in general, since that 
property calls for a wavefunctional which, on large scales, is 
Gaussian in the field strengths rather than the gauge fields. On 
the other hand, also unlike the other proposals, the Gaussian 
wavefunctional is tractable analytically. 

Efforts in this direction were spearheaded by Szczepaniak 
and Swanson [9, 33]. They used a Coulomb gauge ground 
state wavefunctional of the form 



*P[A] =^exp 



d u k 
(27rp 
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The proposal was further developed in ref. [10], where the 
contribution from the Faddeev-Popov determinant was in- 
cluded at one-loop order. The field-independent function G)(k) 
was determined from a gap equation obtained by minimizing 
the energy expectation value. The gap equation depends on 
the so-called ghost dressing function d(k), which is defined in 
terms of the expectation value of the inverse Faddeev-Popov 
operator 4 



d D xe ikx 



r.OM 



= 8 



ab 



k 2 



(39) 



and the Coulomb form factor, f(k), defined by 



/(*) 



$d D xe lkx l?¥\ 



(-D-V) 



!d D xei k X{ ^y_^ )xa . A 



(40) 



In terms of d(k) and f(k) the expectation value of the 
Coulomb kernel in eq. (35), which determines the Coulomb 
potential V, is given by 



V(k) = I d D xj kx (*¥\F ab (x,Q,\A])\V) = 8 



abf(k)d 2 (k) 



k 2 



(41) 

Finally, inclusion of the Faddeev-Popov determinant at one- 
loop order introduces dependence on the function 5 (k = k'/\k\) 



Nc 
2 



d q 
~(2k) 2 



[l-{k-qf] 



d(q)d(q — k) 
(q-k) 2 ■ 



(42) 



which is related to the expectation value of ^ . In ref. [10] 
X(k) (there denoted by F(k)) was derived in context of the 
gap equation, while the explicit representation of i in terms 
of x (k) was derived by Reinhardt and Feuchter in ref. [ 1 2] (cf . 
eq. (47) below). 

The set of coupled Schwinger-Dyson equations for 
X(k),d(k),f(k) and <X>(k) is UV divergent and requires 
renormalization. In the variational approach this is achieved 
by adding relevant and marginal counter-terms to the Hamil- 
tonian and, if needed, renormalizing the functional measure. 
The latter was obtained in [10] and reads 



X(k)^x{k,n)=I x (k)-I x (n) 



(43) 



where I%(k) is given by the right hand side of eq. (42). In 
[10] the renormalization program was, however, not fully im- 
plemented. In particular a Hamiltonian counter-term propor- 
tional to / All, which defines the c\ renormalization constant 



4 As shown by Reinhardt [34], in Coulomb gauge the inverse ghost form 
factor dT 1 (k) has the meaning of the dielectric function of the Yang-Mills 
vacuum, and the horizon condition 

d~ 1 (0)=0 (38) 

therefore implies that the Yang-Mills vacuum is a dual superconductor. 

5 For later use, we present all explicit expressions in D = 2 space dimensions 
and for the color group SU(Nc) [13]. 



(cf. eq. (52) below), was omitted and thus only an approx- 
imate low-energy solution could be obtained. It was found, 
however to be qualitatively consistent with the results of [9] 
that used the J? = 1 (x(k) = 0) approximation. This hints that 
within the one-loop variational approach, contributions from 
the FP operator may be accounted for by the gaussian wave- 
functional itself, with an appropriate choice of the gaussian 
parameter co(k). Such a possibility was rigorously demon- 
strated by Reinhardt and Feuchter [12] (cf. eq. (46) below 
and the discussion that follows). 

Inspired by the wavefunctional of a spinless particle in an 
s-state of a spherical potential Feuchter and Reinhardt in [11] 
suggested to use the ansatz 



VJW] 



exp 



d 2 k 



(2k 



p (o(k)Af(k)A?(-k) 



(44) 

which has a number of technical advantages: The factor of 
t fl\A\ in the integration measure (eq. (36)) cancels against 
J* [A] from the square of the wavefunction and thus drops 
out from the calculation of equal-time vacuum expectation 
values. As a consequence Wick's theorem can be applied di- 
rectly, and in particular (o(k) appearing in eq. (44) is found to 
be directly related to the static gluon propagator 

(Af(fc)A^)) = {2nf8 2 {k + q)8" b ^^. (45) 

In ref. [12] Reinhardt and Feuchter considered a general wave- 
functional of the type 



ValA] 



J/ 



'[A] 



exp 



d 2 k 



^A(-*)(Ha(*)A(*) 



(46) 

In the one loop approximation they showed that the Faddeev- 
Popov determinant, eq. (33), can be represented as 

J^L A ?(-k) X (k)At(k) 

where x(k), thereafter referred to as the curvature, is given by 

5 2 ln J 



c /[A]=exp 



(47) 



8 ab X(k) = -\[d 



2 xe ikx (W a 



8A"(x)8A b {0) 



, (48) 



which, to the order of approximation considered, after renor- 
malization, coincides with the one given in eq. (43). Combin- 
ing eq. (46) and eq. (47) leads to 



^ a [A]-^Kexp 



d 2 k 



A(-jfc) [co a (k)-2ax(k)]A(k) 
(49)" 



and establishes equivalence, at a one-loop level, between the 
ansatz of the Indiana group eq. (37), which corresponds to 
a = 0, and that of the Tiiebingen group eq. (44), correspond- 
ing to a = 1 /2. 6 



6 The value of a does not matter in the one-loop approximation considered 
here. It will, however, become relevant for calculations at higher loop or- 
der. 
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However, using equivalent variational ansatze did not lead 
to the same results for the correlation functions, d(k), f(k), 
%{k), <a(k). This is because the approaches of the Indiana and 
Tubingen groups differ in i) the approximation scheme used 
to evaluate the expectation value of the Hamiltonian and ii) 
the renormalization scheme. While the Tubingen group fully 
includes the Faddeev-Popov determinant to the order consid- 
ered, the Indiana group set Jf — 1 throughout ref. [9] and ne- 
glected J* in the Coulomb term in the numerical calculations 
of ref. [10]. (In the analytic calculation of ref. [10] J* was, 
however, fully included.) Also, while the Indiana group con- 
siders the one-loop corrections to the Coulomb form factor 
f(k), the Tubingen group employs the d(k) = 1 approxima- 
tion in the equation for f(k). 

Ref. [10], in which the renormalization program was not 
fully implemented, missed a Hamiltonian counter-term pro- 
portional to / All, which defines the c\ renormalization con- 
stant (cf. eq. (52) below). The existence of this term was real- 
ized by Feuchter and Reinhardt [11], who carried out the com- 
plete renormalization program. The c\ counter-term missed 
in [10] plays an important role in determining the IR proper- 
ties of the wavefunctional, as realized by Reinhardt and Ep- 
ple [35], and will be crucial for the investigations given in the 
present paper. Therefore throughout this paper we will use the 
fully renormalized approach of the Tubingen group [11, 35] . 

For later convenience we define 



W(k) = co(k)- X (k), 



(50) 



where <X)(k) corresponds to the wave functional in eq. (44), 
and write the wave functional of eq. (44) in the form 

^L A (-k)W(k)A(k) . 



^coAl^^exp 



(2*) 2 



(51) 

The fully renormalized gap equation for co, which ultimately 
determines W, reads [11, 35] 

co 2 (k) =k 2 + x 2 (k)+c 2 + AI {2) (k) + 2 X (Jfc) [A/ ( 1 > (jfc) + c x } , 

(52) 

with 

A/W(Jt)=/W(jfc)-/W(0), 



/«(*) = ^ f ^ 



(2k)- 



(k-q) 2 V(q-k) 



W(q)-W(k) 
(0(q) 



(53) 

and V (k) given by eq. (41). The gap equation, together with 
eq. (43) and the Schwinger-Dyson equations for the ghost 
form factor, 

d-\k)=d-\ l i)-{I d (k)-I d ( i l)), 



ld{k) 



N f 



d q 



2 J (2n) 2 
and Coulomb form factor, 

/(*)= /00 + ('/(*) -W) 



[l-(K-q) 2 } 



d(g-k) 
(0{q) (q-kf 



(54) 



N r 



form a closed set of coupled integral equations for %,d,f and 
CO. In the gap equation (52), c\ and C2 are (finite) renormal- 
ization constants. For the critical solution, where one imposes 
the horizon condition for the ghost dressing function, eq. (38), 
both Co(k) and x{k) are infrared divergent, which implies that 
the transverse gluon propagator vanishes at k — s- 0, while [35] 



W(0) = hm((o(k)-x(k)) = c 1 . 

k^0 



(56) 



So even when enforcing the horizon condition, the quantity 
c\ = ft)(0) is undetermined and may be taken to be either in- 
frared finite or zero. However, a perimeter law of the 't Hooft 
loop requires c\ = and this value is also favoured by the 
variational principle [35]. Furthermore, for ci = 0, in the IR 
limit k — >• 0, the wavefunctional eq. (51) becomes independent 
of the gluon zero mode which agrees with the behavior of the 
exact vacuum wavefunctional in 1 + 1 dimensions [36], and 
corresponds to the so-called ghost loop dominance in higher 
dimensions [37]. But although there is strong evidence to fa- 
vor ci = ft)(0) = 0, our numerical studies in Section IV B will 
also look at the case of a non-zero, but small, value for W(0). 
The renormalization parameter ci, on the other hand, has no 
influence on the IR or UV behavior of the solutions of the gap 
equation (52). Only the mid momentum regime of <X)(k) is 
weakly dependent on C2 [11]. Since we are mainly interested 
in the IR properties we will put co — throughout this paper. 

The set of coupled integral equations can be solved ana- 
lytically in the IR (for the critical solution) using the power 
law ansatze [11, 38] while the full numerical solutions of the 
above equations were given, for D = 3 space dimensions, in 
[11, 39, 40]. For D = 2, the numerical solution was presented 
in ref. [13] and it will be used in Section IV B for comparison 
with lattice simulations. 

One criticism that can be leveled at the Coulomb gauge pro- 
posal is that it is not clear how it could ever lead to an area 
law falloff for spatial Wilson loops. In order to address this 
issue, a modified version of a Gaussian ansatz, which incor- 
porates monopole configurations, has been proposed by Mat- 
evosyan and Szczepaniak [41]. Furthermore, recently [42] 
Campagnari and Reinhardt have developed a method which 
allows to use non-Gaussian wavefunctionals in the variational 
approach. Specifically, a wavefunctional containing vertices 
with up to four gluon fields was considered. Tests of these 
modified versions are, however, deferred to future investiga- 
tions. 



III. THE MEASUREMENT METHOD 



We begin with the identity 



^kwi = zj du j n n^M)-^)] \ e 



-s 



(57) 



where, in the infinite volume limit, is the ground state of 
the operator H, defined via the transfer matrix T = exp [— Ha t \ , 
with a, the lattice spacing in the time direction. In the con- 
tinuous time limit, H is the Hamiltonian of the lattice gauge 
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theory. Now consider a finite set of lattice configurations 
^ = {U^"\x),m = 1,2, ...,M} at a fixed time, and define 

~ M r { 2 1 

z =l ° u nn 5 ^(x,o)-f/ A (m) ( X )] e - s ( 58) 

m=l 7 [ x jfe=l J 

This is the partition function of a statistical system in which 
the lattice configurations at time t = are restricted to the set 
^ . The rescaled wavefunctional 

^[t// B) (x)] 

{fl nLi5[f4( X ,o)-c/f ( X )]} e - S 
~ e^i/w {n x nLi5[^(x,o)-f/, (m) ( X )]} ( 



,-5 



(59) 



has the interpretation as the probability P„ that, in this statis- 
tical system, a lattice configuration on the t = time-slice is 

equal to the n-th configuration U- (x) £ ^ in the given set. 

The probability P n can be computed numerically by a mod- 
ified lattice Monte Carlo simulation. All links at t ^ are 
updated in the usual way, which for the SU(2) gauge group 
with the Wilson action is a simple heat bath. On the t = 
plane, however, one of the M configurations from the set 
is selected at random, and then accepted or rejected by the 
Metropolis algorithm. Let N„ represent the total number of 
times, in a given simulation, that the n-th configuration in the 
set is selected by the Metropolis algorithm, with N lol the total 
number of updates of the t = plane. Then 



N lo ,^°°Ntot 



Nn 



(60) 



Since *Fo[i/(™)] is simply a constant rescaling of 1'o[t/^], it 
follows that the relative amplitudes of the vacuum wavefunc- 
tional 'I'o in the set % are given by 



%[UW] N„ 
— =- — ~ ( tt = lrm — 

*p2[[/(m)] N t0 ,^°°N m 



(61) 



Now suppose we have some theoretical proposal for the 
Yang-Mills vacuum wavefunctional 



theory 



(62) 



If the proposal is correct, i.e. ^theory = ^o, and we make a 
plot of 



log 



Ntn 



vs. R[U { 



(63) 



then the data points should fall on a straight line, with slope 
equal to one. 

The method just described was introduced and applied in 
refs. [15-17]. In that early work, however, the simulations 
were carried out on small lattices and relatively small values 



of j3 = 4/g 2 , while comparison to theory was limited to sim- 
ple wavefunctionals, resembling (6), inspired by the strong- 
coupling expansion. It is now possible for us to greatly im- 
prove on these previous studies. 

In this investigation we will consider sets of three different 
different types of configurations: 

• Abelian plane waves with fixed wavelength L and vari- 
able amplitude 



U[ m) (n l5 n % ) = ^/l-(aW(n 2 )) 2 l2 + ia [m) (» 2 )cJ3 
t/ 2 (m) («i,« 2 ) = 1 2 

a' m '(« 2 ) = -yja + ymcos 



(64) 



where m = 1 , 2, . . . , m max with L the lattice extension and 
a, 7 some constants. 

• Non-abelian constant configurations, variable ampli- 
tude: 7 



U[ m \n u n 2 ) = ^/l-( fl ('»))2i 2 + / fl ('») (Tl 

t/ 2 (m) («!,« 2 ) = V /l-( fl W)2l 2 + /flW<T 2 



» _ 



a + ym 



-,1/4 



20L 2 



(65) 



■ Non-abelian constant configurations, fixed amplitude, 
variable "non-abelianicity" specified by an angle 9,„ 

u[ m '(ni,n 2 ) = y/ 1 -a 2 l 2 + /ac7 1 

f/ 2 (m) (« 1 ,« 2 ) = v / l-a 2 l 2 

+ia (cos (0„,)cti +sin(0 m )(7 2 ) 
m = y(m— (66) 



IV. RESULTS 

Since the measurement method in the previous section re- 
lies on a lattice regularization, we must apply this regulator to 
the vacuum wavefunctionals under study. Let us begin with 
^GO- The proposal is that 



-log[V 2 GO [A]}=R GO [A]+R , 
where Rq = — \og(jY 2 ), and in the continuum 

R G o[A} = -^ Jd 2 xJd 2 yB a (x) 

-, ab 



(67) 



V-z> 2 -V 



B b (y). (68) 



7 The factor of 20 in the definition of a'" 1 ' is an arbitrary scaling of the pa- 
rameters, which could of course be absorbed into a,y. 
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In the special case of abelian plane waves with A"(x) 



A i (x)8" 3 , A 2 (x) = 0, we have the simpler expression 



1 



R G o[A} = ^ d 2 x d 2 y {d 2 A l ) x 



(*Ai), 



(69) 



The engineering dimension of the kernel, in 2+1 spacetime 
dimensions, is also inverse length. We now latticize the theory 
and absorb dimensions into a lattice spacing a, with 



Ai(*)-> -A Ll (x) , d 2 ^ -<?L2 
a a 



d 2 x - 



Si 



pa ' 



(70) 



where di is the lattice finite difference operator, and all of the 
other subscript L quantities are dimensionless. All factors of 
a cancel in R[A], and the result is 



R G o[A] = ^Y^{d L2 A LX ) x 



1 



-V 2 +m 



{d L 2A L l)y 



(71) 



A. The GO and KKN wavefunctionals for abelian plane waves 

Now we specialize to the lattice abelian plane wave con- 
figurations listed in the previous section (lattice sites are x = 
(ni,n 2 )) 



3 u^im^-u^im,^) 



2i 



U 2 v, (n u n 2 ) = l 2 



Ji)/ \ 2 , ; / 2nn 2 

A u ("2) = - V « + YJ cos ( 



I — 



(72) 



Substituting these configurations into R[A], the result is 

R GO [U {]) ] = 2(a + yj)co GO (k 2 ) , (73) 

with 

J3 k 2 



(Oco{k 2 ) 



k- 



? 2 \Jk 2 + m 2 ' 



(74) 



and where k and m are the momentum and the mass parame- 
ters in physical units, i.e. k 2 = k /a 2 ,m = m^/a. 

The same regularization applied to the KKN wavefunc- 
tional yields, for the abelian plane wave configurations, 



Rkkn[U U) ] = 2(a + yj)(o K K N {k 2 ) 



(75) 



with 



G>KKN{k ) = 



k 2 



(76) 



The theoretical values for (o{k 2 ) are to be compared against 
the data obtained from the numerical simulation. For a given 
lattice coupling fa of the Wilson action, at a given lattice size 
L corresponding to a value of k given in eq. (72), we obtain 
from the numerical simulation described in the previous sec- 
tion the values 



Then a>Mc(k 2 ) is obtained from a best linear fit of 



2{a + yn)(o M c{k 2 



(77) 



(78) 



to the data points {r„}. Figure 1 shows a typical plot of r„ vs. 
2(a + yn) at fi E = 9 and L = 24; (OMc{k 2 ) is the slope of the 
line (best linear fit) shown. The values for a, y used at each 
Pe and L are listed in Table III of the Appendix. 

The theoretical expressions for ©go(^ 2 ) an d ®kkn {k 2 ) in- 
volve two dimensionful parameters m and g 2 . Once these pa- 
rameters are chosen, the results can be compared with the data 
obtained for 0)mc(^ 2 ) on an Y lattice, providing the dimension- 
less squared momentum k on the lattice is converted into 
physical units k 2 = k ja 2 using the lattice spacing a. For a 
choice of lattice coupling Pe, the lattice spacing in physical 
units is given by 



(79) 



where O/, = Ol{Pe) is the D = 3 dimensional string tension 
in lattice units, and a is the string tension in physical units. 
On grounds of tradition, we make the arbitrary choice a = 
(440 MeV) 2 . 

Figure 2 is a plot of a)Mc{k 2 ), extracted from a best fit of 
the data to eq. (78). Each data point is obtained at a particular 
Pe =6,9 or 12 on a given lattice of extension L, with L = 
16,24,32,40 or 48, and the wavelength of the plane wave on 
each lattice is the largest wavelength A = L available. This 
plot also displays the two theoretical curves 



(Oco{k 2 ) 

(OKKN{k 2 ) 



1 



g 2 Vk 2 -^' 
1 



-mr 

a 



(80) 



with the parameters g 2 and m obtained, for each curve, from 
a best fit to the data points. Observe that in this range of mo- 
mentum, the difference between the two fitting functions is 
essentially negligible, and in fact only becomes noticeable for 
k 2 > 4 GeV 2 . 
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abelian plane wave, 




lattice spacing at each [$e- The GO result is within 5% of that 
value, and the KKN result is almost exactly right. 



FIG. 1. A typical plot of the data for -log(N n /N tot ) at jS £ = 9 and 
lattice extension L = 24, vs. the factor 2(a + yn) associated with the 
amplitude of the ;i-th configuration. The straight line is a best linear 
fit, and the quantity 0) MC (k 2 ) is the slope of that line. 



a. 

o 

3 




p 2 (GeV 2 ) 



FIG. 2. Cumulative data for (Omc vs - P 2 m physical units, on lattices 
of extensions L = 16,24,32,40,48, and Euclidean lattice couplings 
Pe = 6,9, 12. The curves labeled "GO fit" and "KKN fit" (there are 
actually two curves, difficult to distinguish from one another), are 
the theoretical values for ft>co(p 2 )> an d ®kkn{p~)> using the param- 
eters of m and g 2 in Table I. The line labeled "Coulomb gauge" is 
obtained from the ansatz for the Coulomb gauge vacuum wavefunc- 
tional fee [A] (eq. 51) as described in Section IV B. 



With the parameters obtained from the fit, we can use di- 
mensional reduction (naively, in the KKN case, as explained 
in section II C) to compute the string tension, and compare it 
with our input value of (440 MeV) 2 . Dimensional reduction 

gives 



(81) 



The parameters g 2 ,m from the best fit, and \fa from obtained 
dimensional reduction, in the GO and KKN cases are shown 
in Table I. The values of y/a should be compared with the 
given value of \fa = 0.44 GeV, which was used to set the 




variant 


m 


2 

8' 


\fd from 
diml red. 


GO 
KKN 


0.771 
0.420 


1.465 
1.237 


0.460 
0.441 



TABLE I. The parameters m,g 2 for the GO and KKN wavefunction- 
als, determined from a best fit to the abelian plane wave data in Fig. 
2, with \fa derived from dimensional reduction. All values are in 
units of GeV. 

The product of m and g 2 , in either the GO or KKN ap- 
proach, determines the string tension a in either approach. 
The dimensionless ratio g 2 /m is an output of the KKN ap- 
proach, where it is predicted to be n. If m and g 2 are de- 
termined from a best fit to the data, then the actual ratio is 
g 2 /m = 2.95. It is not clear, at this stage, whether this small 
discrepancy is significant, or should just be attributed to devi- 
ations from the continuum scaling due to a finite lattice spac- 
ing. 



B. Tests of the Coulomb gauge wavefunctional 

To test the wavefunctional eq. (51), we first have to trans- 
fer it to the lattice. We begin by rescaling the gauge field 
Aj i ^ A(/g so that a prefactor g~ 2 appears in the exponent 
of eq. (51), and A, (x) has engineering dimension of a mass. 
With these conventions, the Fourier transformed kernel (o(k) 
and curvature %(k) also have dimensions of mass. 

Next we latticize as in eq. (70) and rescale the gauge field 
again to obtain the dimensionless field 8 A k (x) = aA c k (ax). For 
Coulomb gauge fixed connections, it is, in principle, impor- 
tant to use the so-called midpoint rule when extracting the 
gauge fields from the lattice links t/j: 



U k {x)=a° k (x)l + ia L k (x)a c 
Al(x + i/2) = -2al(^-n(4(x)). 



(82) 



As compared to simpler prescriptions such as eq. (72), we 
have two modifications: 

1. The shift in the argument on the lhs ensures that the 
resulting lattice connection is exactly lattice transversal 
if the link fields are, 



V-A(i)=£ Aj{x + J)-Aj(x) 



0. 



After Fourier transformation, the shift leads to a phase 
factor in the connection which affects general observ- 
ables but happens to drop out in the (quadratic) expo- 
nent/? [A] tested here. 



' Throughout this section, we will denote dimensionless lattice objects with 
a caret. 



2. The 77-correction in eq. (82) comes from the SI/ (2) al- 
gebra for parallel transporters over a finite distance a, 

, . arccos t , 
VI -f 2 

It is only relevant for very strong fields far from the con- 
tinuum limit. (In our numerical studies, the correction 
never exceeded 5%.) 



After Fourier transformation 



A?(*)=J>-*% e (*) 



(83) 




where k, — (2n/L)£j (with — L/2 < £j < L/2), a simple calcu- 
lation leads to the lattice version of the CG wavefunctional, 



2 3 



= lc=l 



W(k) 



[(o(k) — x(k)] ■ 



Ra 



(84) 



Notice that the dimensionless momentum argument in the nu- 
merical continuum solution of the gap equation is k/g 2 , so 
that its lattice counterpart becomes 



-2 (% \ 
h = — ~ sm 



(85) 



To complete the lattice transcription, we only have to find an 
expression for the function 



h(fi) = a(P)i 



(86) 



where j3 =A/(agfy is the usual lattice coupling for SU (2) MC 
simulations in D = 2 + 1 . From high precision measurements 
of the string tension in D = 2 + 1 [43], the best fit in the scaling 
window J3 G [3,12] is 



(7a 



b (< c 

Fl 1 + ? 



with coefficients b « 1.788 and c « 1.414. From this, 



2 16 
a = a a = a 



16(7 



si 

From the leading terms of order 
and therefore 



), we find b = I6a/g 4 



h(fi) = ag 2 = ^ = VW)^= 



4 / c" 



\[a ^b 

c= 1.414. (87) 



This completes the lattice transformation of the Coulomb 
gauge wavefunctional. 

Let us first look at the non-Abelian constant configurations 
(65). The corresponding lattice connection has the special 
colour structure A- ~ 8[, but is otherwise constant in space, 



FIG. 3. The exponent R from the variational approach eq. (88) plot- 
ted against the lattice data for — In 1 ? 2 for one set of non-Abelian 
constant configurations, choosing ffl(O) = c\ as fitting parameter 
(d =0.1165). 



i.e. Fourier transformation projects out the zero frequency 
contribution, 

££k(*)i 2 ~s* i0 . 

;=ic=i 

The final result for the exponent in the wavefunctional 

non-Abelian constant con- 
figurations, 



Rcg[U [ ' 



8L 2 arccos 2 



Rn 



8L Z {a [m> ) 2 -W(0)+R 



(88) 

where the approximation in the second line comes from dis- 
carding the 77-correction in eq. (82). 

From eq. (56), the quantity ft)(0) is given by the (finite) 
renormalization constant c\ and, as already mentioned in 
sect. II D, the energetically preferred value is c\ = 0, which 
is also required for a perimeter law in the 't Hooft loop 
[35]. Obviously, with this choice of renormalization constant 
(fl(0) = c\ = the Coulomb gauge wavefunctional cannot ac- 
count for the constant non-Abelian gauge field configurations. 
Whether this failure is important remains to be seen. At least 
it does not necessarily imply that the Coulomb gauge wave- 
functional is a bad approximation to the true vacuum wave- 
functional since constant configurations form a set of measure 
zero in field space. One could give up the preferred value 
C\ = and choose ft)(0) = c\ as a fitting parameter, cf. fig. 3. 
This gives reasonable agreement with the lattice data for one 
set of constant non-Abelian configurations but does not cure 
the general problem. From the results presented in Sec. IV C 
below, it will become clear that constant non-Abelian gauge 
fields can only be accounted for if we include quartic terms 
~(4x A) 2 in the exponent of the wavefunctional. The use of 
such non-Gaussian wavefunctionals in the variational princi- 
ple has recently become feasible [42], but the solution for the 
wavefunctional has not yet been determined explicitly up to 
quartic terms in the exponent. 



For these reasons, we will use the energetically favored 
value ftJ(O) = c\ = in the following. We will now show 
that the Coulomb gauge wavefunctional does a good job for 
Abelian plane waves of the type eq. (64). In this case we have 
carried out simulations at j3 = 6 on a fixed lattice volume of 
extension L = 24, and varied the amplitude of the plane waves, 
at given wavelength L/M, according to 



U[ m \n u n 2 ) = v /l-( fl ('")(« 2 )) 2 l 2 + /fl(" , )(« 2 )(J3 

f/ 2 (m) (m,n 2 ) = l 2 

(m), \ 1 / (2Kn 2 M\ 

a 1 '(«2) = -s/m k m cos I — - — J , 



(89) 



where m = l,...,m max , with k m = 1.4,0.45,0.17,0.09,0.036 
at M = 1,2,4,8, 12 respectively. The connection is Abelian, 
A- ~ 5 t3 , with a harmonic spacetime dependence in the y- 
direction; the corresponding wavenumber is proportional to 
the parameter M in eq. (89). After Fourier transformation the 
general result (84) takes a fairly complicated form 



L/2 

R CG [U {m) }=R + 4 £ W(p n ] 

«=-L/2+l 



xsgna w (r) • arccos J I - (d m > ) 2 {r 



Pn 



HP) 



sin 



% 

— n 
L 



(90) 

This can be simplified considerably, if the 77-correction in the 
definition of the connection, eq. (82), is discarded. Then the 
sums in eq. (90) can be performed explicitly and we obtain a 
much simpler expression 



Rcg[U { 



R + 2cM-mK M -Co(p M ) 



(91) 



where cm = 2 for the highest frequency M = L/2 and cm = 1 
otherwise for L even (L = 24 in this case). From eq. (91), it 
is obvious that the plane wave configuration tests the kernel 
W = (0 / g 2 — % / g 2 exactly at the lattice momentum ~p M which 
corresponds to the frequency of the plane wave. 

Figure 4 shows the result of the numerical evaluation of 
eqs. (90), (91) against the lattice MC data for Abelian plane 
wave configurations of varying wavenumber and amplitude. 
As can be clearly seen, the individual plane waves with fixed 
wavenumbers M and varying amplitude fall on a straight line, 
but the slope of that line differs from unity. (We have chosen 
the solution W(k) of the variational problem with the preferred 
renormalisation constant c\ = 0.) Morever, the slopes of the 
lines vary slightly with M, i.e. effectively with the momen- 
tum picked by the plane wave: For the smallest momentum 
M = 1, we find a slope of 1.19, which decreases down to 1.02 
for M = 2, and then increases again up to 1 .52 for the largest 
momentum M = 12 representable on a L = 24 lattice. If we 
relax the condition on the renormalisation constant c\ and take 
it as a free parameter, we observe that the spread in the slope 
between the various wave numbers is increased, which is an- 
other hint that the choice c\ = should be preferred. 




FIG. 4. The exponent Rcc from the variational approach eq. (90) 
plotted against the lattice data for — In I* 2 for the plane wave config- 
urations with wavenumber M £ {1,2,4,8, 12}. The lattice data was 
taken with lattice extension L = 24 at j8 = 6.0. 



Since the plane waves test the kernel (o(k) at varying mo- 
menta, we can use a fit to the MC data as explained in the 
previous section to find a numerical estimate a>Mc(k)- In the 
Coulomb gauge wavefunctional, this quantity corresponds to 
~o5(k) = g~ 2 (co(k) — X{k))- After rescaling to physical units 
(see eq. (86) and below), the result is plotted along with the 
values obtained by numerical simulation, G>Mc{k), in fig. 2. It 
is evident that the variational solution for W(k) fits the MC 
data very well, at least in the infrared region for momenta up 
to k w 1.3GeV. For larger momenta, W(k) starts to deviate 
and becomes slightly larger than the numerical estimate, but 
at most by a few percent within the phenomenologically rel- 
evant mid-momentum regime. (For very large momenta not 
plotted here, c5(k) ~ k is exact by asymptotic freedom.) 



C. Non-abelian constant configurations: fixed amplitude, 
variable "non-abelianicity" 

For general non-abelian configurations we have, in a lattice 
regularization, 



RgoP 



where 



(»)l 



-Lin 




ah 



1 



-D 2 -Ao+m? 



B"{x) = -Tr[U{P x )o")] 



B b (y) 



(92) 



(93) 



with U(P X ) a product of links around a plaquette, starting with 
a link at site x. The lattice covariant Laplacian, in the adjoint 
representation, is given by 



k=\ 

■ ' ' O a U k {x)o b ul{x) 



(94) 
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In terms of the parameters g 2 7 m in the GO row of Table I, we 
use j3 = 4/ (g 2 a) and = ma, where a is the lattice spacing. 
For comparison with the Monte Carlo data generated at the 
lattice coupling fa of the Wilson action, we determine a from 
eq. (79). It is important to note that while we expect j3 /fa — > 1 
in the continuum limit, this ratio need not be exactly equal to 
one at any finite fa- 
in the same way, the latticized "hybrid" wavefunctional is 



B b (y) 



(95) 

with p , ml determined using the parameters g 2 , m in the KKN 
row of Table I, and the lattice spacing from eq. (79). 

We will consider first the configurations of eq. (66), with 
fixed amplitude and variable "non-abelianicity" specified by 
the 9 parameter. If the amplitude is chosen small enough, 
then — D 2 — Ao is negligible compared to m 2 , and the kernel 
reduces to 




1 



\^-D 2 -Xo+m 2 
for the GO wavefunctional, and 



= — S XY S 
m 



ab 



(96) 



ab 



2m 



ab 



(97) 



for the hybrid. This is the dimensional reduction limit, and in 
either case, for the configurations (66), R[U] « (A[ X A2) 2 , or 



R 



GO.hybrid 



[C/M]c<sin 2 (e„) 



(98) 



For the Coulomb gauge wavefunctional, however, /?[{/]« 
A 2 +A\, and hence, since the amplitudes of A\ and A2 are 



fixed in the set (66), 



R CG [U {n) ]~W{0) 



(99) 



independent of the angle 0„. If ft)(0) = 0, which seems opti- 
mal for agreement with the plane wave data, then Rcg would 
also be independent of the amplitude of the gauge fields. 
However, it is important to recall that the Coulomb gauge 
wavefunctional should not be evaluated outside the first Gri- 
bov horizon. So even if ft)(0) = 0, the restriction to the Gri- 
bov region amounts to a cutoff in the amplitude of non-abelian 
constant configurations. 

The Monte Carlo simulation was carried out on a 12 3 lattice 
at fa — 6, with the t = configurations chosen from 



Faddeev-Popov operator, we have checked that these lattice 
configurations are all inside the first Gribov horizon. 

In Fig. 5 it can be seen that the logarithm of the wavefunc- 
tional is indeed proportional to sin 2 (0), as one would expect 
from the GO and hybrid wavefunctionals in the dimensional 
reduction limit. The data does not seem to be compatible, 
however, with the ©-independence (99) of the CG wavefunc- 
tional (51). 

We recall that if ^[U] = exp[—jR(U)] is the true vacuum 
state, then the data points for — log(iV n /iVr) vs. R[U n ] should 
fall on a straight line, with unit slope. Plotting the data for 
— log(N n /NT) against Rqq [U n ], as in Fig. 6, we find the slope 
obtained from a linear fit through the data is indeed close to 
unity. In the GO case the slope is 1.02(6); a similar analysis 
for the hybrid wavefunctional results in a slope of 1.12(7). 
Some numerical details concerning the simulations are found 
in the Appendix. 
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FIG. 5. Dependence of —log(N n /Nj) on the "non-abelianicity" of 
the non-abelian constant configurations, determined by sin(6„). 
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Uf> = Vl-« 2 l 2 + /a(cos(0„)c7i+sin(0„)cr2) (100) 

with a = 0.193, and 9„ = (n - 1)^/38. By explicitly calcu- 
lating numerically the lowest lying eigenvalues of the lattice 



FIG. 6. Plot of —\og(N„/NT) vs. Rqq f° r the non-abelian constant 
configurations with variable non-abelianicity. The straight line fit has 
slope = 1.02. 
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D. Non-abelian constant configurations: variable amplitude, 
maximal "non-abelianicity" 

We now consider the non-abelian constant configurations 
of maximal "non-abelianicity,", i.e. = 7C/2, which are the 
configurations of eq. (65), with index m running from 1 to 20. 
All Monte Carlo calculations were carried out on lattices of 
volume 32 3 at Pe = 6,9, 12, and the corresponding values of 
P,niL at each Pe are given in Table II, where the values for 
the hybrid wavefunctional are taken to be the KKN values, 
since the hybrid reduces to the KKN form on abelian config- 
urations. The test of the GO and hybrid wavefunctionals is 
to see whether or not the data points for — log[N„/N tot ], when 
plotted against R\U^ n '], fall on a straight line whose slope is 
close to unity. 



& 


/3(GO) 


m L (GO) 


j8 (KKN) 


m L (KKN) 


6 


4.73 


0.445 


5.60 


0.242 


9 


7.43 


0.283 


8.80 


0.154 


12 


10.19 


0.207 


12.07 


0.113 



TABLE II. Values of /3 , m L for the GO and KKN wavefunctionals at 
each Pe, derived from the g 2 ,m parameters in Table I and the lattice 
spacings a, at (5e = 6,9, 12. 
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FIG. 7. Plot of —log(N n /Nj) vs. Rqq f° r non-abelian constant con- 
figurations, maximal non-abeliancity, at Pe = 6, L = 32, a = 2, y = 
0. 15 In this case the straight line fit has a slope = 0.98. 

An example of the — log[N n /N tot ] vs. /? G0 [!/(")] data at 
Pe = 6 is shown in Fig. 7, for the choice a = 2, y = 0.15. 
Although the data is nicely fit by a straight line which has a 
slope close to unity, this fact must be interpreted with cau- 
tion because, since the number N„ falls off exponentially with 
Rqo[U(">], the range of R must necessarily be kept small; typ- 
ically AR ps 4 — 5. This could mean that the tendency of the 
data to lie on a straight line is misleading, and perhaps we are 
simply looking at the tangent of a curve. It is therefore neces- 
sary to extract the slope of the straight line over small intervals 
centered around points over a wide range of R. The question 
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FIG. 8. Slopes for the GO wavefunctional vs. R, at Pe = 6,9, 12 and 
L = 32, using the values of g 2 , m derived from the abelian plane wave 
fit. 
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FIG. 9. fe=12 calculation, for both types of wavefunctionals. 



is whether those slopes are constant, in which case the linear- 
ity hypothesis is verified, or whether they vary significantly as 
R increases. This is the motivation to calculate — \og\N n /N to t\ 
in sets of twenty configurations, using different values of the 
parameters (a, 7) in each set. The parameters we have used 
are shown in Table IV of the Appendix. 

Figure 8 is a plot of the slope vs. R at Pe = 6,9, 12, where 
the value of R at each data point is the midpoint of the range in 
which the slope was computed. Things are not perfect; there 
is some slight variation in the slope with R, there is a little 
variation with j3, and the values of the slope are not exactly 
one (they seem to be closer to 1.1 at the large R values). On 
the other hand, we have made no claim that the GO wavefunc- 
tional is exact, nor is asymptotic scaling exact at these lattice 
couplings. The point is that scaling is not bad, and the slopes 
are fairly close to unity over a large range of R, using g 2 ,m 
values that were extracted from fits to a completely different 
type of lattice configuration (i.e. abelian plane waves). 

Results for the hybrid wavefunctional turn out to be quite 
close to those of the GO wavefunctional. The values for Pe = 
12, for both types of wavefunctionals, are shown in Fig. 9, 
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with similar agreement at the two other Pe values. 



E. The ghost propagator and the Coulomb potential 

Because of the equality (4) of the vacuum wavefunctionals 
in temporal and Coulomb gauges, when evaluated on trans- 
verse (V • A = 0) gauge fields, equal-time expectation values 
in Coulomb gauge can be derived from 



P E =9, L=32 



(Q)= DAQ[A]8(V-A 



[A]%[A] 



(101) 



and we may use for either of the temporal gauge proposals, 
I'gOj ^hybrid, or the Coulomb gauge proposal ^cg to calcu- 
late such objects as the ghost propagator 



G(R) 



1 



\V V-D[A] /X 
and the color Coulomb potential 9 

V C (R) = - 



(102) 



lx-y=R 




x-v\=R 



In eq. (101) there is an implicit restriction of the integration 
domain to the Gribov region. For the Coulomb gauge wave- 
functional *¥cg [A] the ghost propagator and the Coulomb po- 
tential are presented in [13]. 

In an ordinary Monte Carlo (MC) simulation, Coulomb 
gauge expectation values are obtained by first generating 
lattice configurations with the usual probability distribution 
exp[— S)/Z, where S is the standard lattice action, transform- 
ing those configurations to Coulomb gauge, and evaluating 
the observable Q in the ensemble of transformed configu- 
rations. In principle the same strategy applies to evaluating 
the right hand side of (101) numerically; the problem in 
that case is to generate configurations with the probability 
distribution ^[I/], and this problem was solved, for the "Pgo 
proposal, in ref. [8]. The simulation method developed in 
[8] is also applicable (although it has not been applied until 
now) to the hybrid proposal. The lattice ghost propagator 
and Coulomb potential were calculated numerically from 
^GO' and compared to the corresponding results in ordinary 
lattice Monte Carlo, in ref. [28]. In that work, however, the 
authors chose j3 = j3g and mi = 4j8(7l/3. In the present 
article the philosophy has changed somewhat. We have two 
parameters with dimensions of mass, g 2 and m, and a scale 
set (arbitrarily) by taking yfo = 440 MeV. Then g 2 ,m are 
chosen to give a best fit to the abelian plane wave data in 
Fig. 2. To compare wavefunctional results with standard 
Monte Carlo res ults w e determine the lattice spacing a, at 
each Pe, from yj <Jl/o, and then p = 4/ (g 2 a) and ml = ma 



' More precisely, for color charges in some representation r, the Coulombic 
potential energy is obtained by multiplying V C (R) by the quadratic Casimir 
C,-, and dividing by the dimension of the adjoint representation. 



rr 
o 




FIG. 10. The ghost propagator derived from standard Monte Carlo 
(MC) simulation at j8g = 9, and the same quantity calculated by sim- 
ulation of the GO and hybrid wavefunctionals, by the technique de- 
scribed in ref. [8]. 



(103) are (h e corresponding dimensionless parameters to use in 



the latticized wavefunctional ^go or 



hybrid ■ 



With the new 



procedure we have ^ j3s, and the obvious question is 
whether this fact will tend to destroy the agreement that 
was found previously, in [28], between ghost propagators 
and Coulomb potentials derived from simulation of ^qq, 
and the corresponding quantities found in ordinary lattice 
Monte Carlo simulations. We would also like to calculate the 
Coulomb gauge ghost propagator and Coulomb potential for 
the hybrid wavefunctional proposal. 

Figure 10 shows the equal-times ghost propagator G(R) 
computed in a standard Monte Carlo simulation on a 32 3 lat- 
tice at Pe = 9. On the same plot we see the corresponding 
results obtained by generating lattices with probability distri- 
bution and *PLj n y by the methods of [8], transforming 
to Coulomb gauge, and evaluating the ghost propagator, in 
each case using the appropriate values of /3 , niL corresponding 
to f$E = 9. It can be seen that the agreement between Monte 
Carlo, GO, and hybrid results is almost perfect. 

The agreement for the Coulomb potential V C (R) is not as 
good. In Fig. 1 1 we display the data from MC, GO, and hy- 
brid simulations, again at [5e = 9, with a cut in the data, dis- 
carding configurations with |V(0)| greater than some bound 
equal to 5, 10,50,300. If we restrict the data set to configura- 
tions with |V(0)| < 5, then the agreement between MC, GO, 
and hybrid results is again almost perfect. Roughly half of all 
configurations meet this criterion. The agreement is still fairly 
good for |V(0)| < 10, which accounts for about 80% of all 
configurations. However, as the cut is gradually removed, the 
Coulomb potential derived from GO and hybrid simulations, 
while roughly linear in R, deviates quantitatively from the MC 
result. But how can there be such a noticeable deviation when 
the ghost propagators agree so accurately, without any cuts at 
all? The explanation probably has to do with a discrepancy in 
the tail of the probability distribution. If two probability dis- 
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FIG. 11. Data for the Coulomb potential at /3g = 9 and L = 32, derived from MC, GO and hybrid simulations, with a cut on the data, discarding 
configurations for which \Vq\ is greater than 5, 10, 50, and 300, respectively. 



tributions agree in their lower moments, but disagree in higher 
moments, then it means that the two distributions agree pretty 
well where the probability is substantial, but disagree in the 
tail of the distributions. That is what seems to be going on 
here. 

What was found already in ref. [28] is that the Coulomb 
potential is quite sensitive to a comparatively small number of 
"exceptional" configurations, in which the lowest eigenvalue 
of the Faddeev-Popov operator — V D is far below the average 
value for the lowest eigenvalue. The reason that such excep- 
tional configurations are relevant for the Coulomb potential, 
but not the ghost propagator, is presumably because the ghost 
propagator involves only one factor of the inverse F-P oper- 
ator, while the Coulomb potential involves two factors. Be- 
cause the inverse F-P operator becomes singular as the lowest 
eigenvalue Xq approaches zero, higher powers of the inverse 
F-P operator (such as the Coulomb potential) will be more 
sensitive to infrequent configurations with exceptionally low 
values of Ao than lower powers (such as the ghost propagator). 
The probability distribution of infrequent configurations is, of 
course, governed by the tail of the probability distribution. So 
our interpretation of the ghost and Coulomb propagatorresults 
is that ^Vqq and l P^ vfc ,. (£/ agree quite closely with each other, 
and with the probability distribution of the true Yang-Mills 
vacuum wavefunctional ^q, in the "bulk" of the distribution. 
The Coulomb potential data suggests, however, there is some 
small disagreement in the tail of the distribution. 

In general, our results for the Coulomb gauge ghost prop- 



agator and Coulomb potential with the new fitting procedure 
for j3 , m agree quite closely with our previous results (based 
on setting /3 = Pe) reported in ref. [28] (for a quantitative 
comparison, cf. [44]). The GO and hybrid results are, once 
again, virtually indistinguishable. Since both choices of pa- 
rameters, and the GO and hybrid wavefunctionals, have about 
the same dimensional reduction limit, our results suggest that 
the quantities we have computed, at the couplings we have 
employed, are mainly sensitive to that limit. 



V. CONCLUSIONS 

We have compared several suggestions for the Yang-Mills 
vacuum wavefunctional to the true Yang-Mills vacuum wave- 
functional in 2+1 dimensions, whose exact form is unknown, 
but whose relative magnitudes in any set of lattice config- 
urations can be obtained numerically. Three types of lat- 
tice configurations were studied: abelian plane wave config- 
urations, non-abelian constant configurations of fixed ampli- 
tude but varying "non-abelianicity," and non-abelian constant 
configurations of maximal abelianicity but of differing wave- 
lengths and varying amplitudes. For purposes of comparison, 
the physical scale was set by taking the string tension to be 
^ = 440 MeV. 

For abelian plane waves, up to the shortest wavelength 
corresponding to p 1 = 2.5 GeV 2 that we have investigated, 
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the GO and Karabali-Kim-Nair proposals are almost indistin- 
guishable, and both agree very well with the values obtained 
for the true vacuum wavefunctional, evaluated on these con- 
figurations. The Coulomb gauge wavefunctional can also fit 
the plane wave data with an appropriate choice of parame- 
ters, providing in particular that the renormalization constant 
c\ in eq. (41) is set equal to zero. Both the GO and KKN 
wavefunctionals reduce to the dimensional reduction form 
exp[— ji J F 2 ] at long wavelengths, and it seems likely that this 
is also true for the Coulomb gauge proposal, in this special 
case of abelian configurations, for the choice of renormaliza- 
tion constant c\ = 0. 

For non-abelian configurations, we have suggested a gauge- 
invariant wavefunctional which reduces to the KKN pro- 
posal for abelian configurations, and incorporates the covari- 
ant Laplacian and eigenvalue subtraction of the GO proposal, 
which we have termed the "hybrid" wavefunctional. Both 
the GO and hybrid wavefunctionals have the dimensional re- 
duction form when restricted to configurations which, when 
expanded in eigenstates of the covariant Laplacian, contain 
only low-lying eigenmodes. Once again, the GO and hy- 
brid wavefunctionals are almost indistinguishable when evalu- 
ated on non-abelian constant configurations, and this is prob- 
ably because they have almost the same dimensional reduc- 
tion limit. We find that the GO and hybrid wavefunctionals 
are in good agreement with the true vacuum wavefunctional 
for non-abelian constant configurations, as well as for abelian 
plane waves. The Coulomb gauge wavefunctional, however, 
which does not have the dimensional reduction property for 
non-abelian lattices, does not seem compatible with the data 
for non-abelian constant configurations, particularly the data 
with variable non-abelianicity. 

The Coulomb gauge wavefunctional has been used to 
compute Coulomb gauge ghost and gluon propagators, with 
results in 2+1 dimensions, reported in [13], indicating a 
Coulomb potential rising almost (but not quite) linearly. We 
have also computed these quantities by direct simulation of 
the GO and hybrid wavefunctionals. The GO and hybrid re- 
sults agree with one another, and almost perfectly with the 
lattice Monte Carlo results for the ghost propagator. The GO 
and hybrid wavefunctionals also lead to an apparently linear 
Coulomb potential and agree very closely with each other. On 
the other hand there is some difference in the GO and hy- 
brid Coulomb potentials in comparison to the lattice Monte 
Carlo results, and this can be attributed to a difference asso- 
ciated with exceptional configurations with unusually small 
values of the lowest Faddeev-Popov eigenvalue. Thus the GO 
and hybrid wavefunctionals would seem to agree with the true 
Yang-Mills vacuum wavefunctional for the bulk of the prob- 
ability distribution, but there would appear to be a small dis- 
agreement in the tail of the distribution. 

The main effort in this article has been to calculate the rel- 
ative magnitudes of the true vacuum wavefunctional on par- 
ticular sets of lattice configurations; namely, abelian plane 
waves and non-abelian constant configurations, and to com- 
pare those results with a number of proposals for the vacuum 
state. We have found that the lattice data for the abelian plane 
waves have been nicely reproduced by all proposals consid- 



ered, while good agreement with the data for non-abelian con- 
stant configurations appears to require wavefunctionals with 
the property of dimensional reduction. 
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Appendix: Numerical details 

Evaluation of Rqq [U ] involves dealing with a kernel 

\ ab 



K" b = 



1 



v/-D 2 -Ao + m 2 



(A.l) 



which, on a lattice of extension L, calls for inverting the square 
root of a 3L 2 x 3L 2 matrix. The numerical evaluation in this 
case can be accelerated using the Zolotarev approximation, 
for which 



1 fl 2 A3 «4 

— = ps a\ H 

y/X X + b 2 l X + b 3 l X + b 4 l 



(A.2) 



where X is a matrix, and the coefficients are given by [45] 

al =0.3904603901 
a2 = 0.0511093775 
«3 = 0.1408286237 
aA = 0.5964845033 
b2 = 0.0012779193 
b3 = 0.0286165446 
b4 = 0.4105999719. 



In fact, what one really wants is the vector 

u a x =K%F» 2 (y), 



(A3) 



(A.4) 



and we found it convenient to compute this vector numerically 
using the Matlab software package. In Matlab, computation 
of the vector u = M~ 1 w, given the matrix M, requires only a 
single line of code: u = M\w. One first defines X = —D 2 — 
Aol + m 2 l to be a sparse matrix, and then sets Y% = X + bit 
etc. The vector u with components u" is then obtained by the 
Matlab statement 

u = ai * 1 + a 2 * (Y 2 \F) + a 3 * (Y 3 \F) +a 4 * (Y 4 \F) , 

(A.5) 
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fe I L = 16 I L = 24 I L = 32 I L = 40 I L = 48 
~6 (005) (OL0) (20,1.5) (30,2.5) (60,3-5) 
9 (3,0.25) (5,0.5) (50,0.7) (10,1.3) (20,1.8) 
12 (2,0.17) (7,0.28) (12,0.53) (20,0.75) (30,1.0) 



TABLE III. Values of a, y used in eq. (72) to generate abelian plane wave configurations with wavelength A = L equal to the lattice extension, 
and fe = 6,9, 12. 



fe l {(«, y)l 

6 (2,0.15) , (15, 0.20) , (32,0.20) , (60,0.22) , (86,0.24) , (107, 0.26) 

9 (2,0.09) , (10, 0.10) , (25,0.13) , (50,0.14) 

12 (1.3,0.06) , (4, 0.06) , (10,0.065) , (20,0.08) , (27,0.083) , (35,0.083) 



TABLE IV. Values of Cf, y used in eq. (65) to generate non-abelian constant configurations with maximal non-abelianicity, on a 32 2 lattice and 
J3 £ =6,9, 12. 



and we finally take the inner product 

R=^F? 2 (x)u a x , (A.6) 

with an implicit summation over lattice sites x and color in- 
dices a. All the matrix operations, including the determina- 
tion of Ao, can be carried out numerically using sparse matrix 
techniques, which results in a considerable savings in com- 
putation time, often by an order of magnitude or more in our 
calculations. We have checked the accuracy of the Zolotarev 
approximation by evaluating R numerically, in several cases, 
without this approximation, and have found the results with 
and without the approximation to differ only at the third sig- 
nificant digit. This is sufficient for our purposes. In the case 
of Rkybrid the formula (A.2) is not directly applicable, and the 
numerical evaluation was carried out without the help of the 
Zolotarev approximation. 

In the Monte Carlo simulations, we set up eight runs each 
time with the same parameters, but different seeds for the ran- 
dom number generator. Each run is itself a number of inde- 
pendent jobs, which we refer to as "cycles", whose results for 
— log(iV n /iVr) are averaged together at the end of the run. At 
the beginning of each cycle the links are all set to the iden- 
tity matrix, except for the spacelike links on the t = plane, 
which are set to the first (n = 1) configuration out of the set of 

(n) 

{W (x,t = 0)} of non-abelian constant configurations. The 
lattice at t ^ then thermalizes for 5000 sweeps with the n = 1 
configuration at t = held fixed. All timelike links are fixed to 
the unit matrix, except for the timelike links at t =L/2, which 
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